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We study interactions between shocks and standing- wave patterns in vertically oscillated layers of 
granular media using three-dimensional, time-dependent numerical solutions of continuum equations 
to Navier-Stokes order. We simulate a layer of grains atop a plate that oscillates sinusoidally in 
the direction of gravity. Standing waves form stripe patterns when the accelerational amplitude 
of the plate's oscillation exceeds a critical value. Shocks also form with each collision between the 
layer and the plate; we show that pressure gradients formed by these shocks cause the flow to 
reverse direction within the layer. This reversal leads to an oscillatory state of the pattern that 
is subharmonic with respect to the plate's oscillation. Finally, we study the relationship between 
shocks and patterns in layers oscillated at various frequencies and show that the pattern wavelength 
increases monotonically as the shock strength increases. 

PACS numbers: 45.70.-n, 45.70.Qj, 47.40.Nm, 47.54.-r 



I. INTRODUCTION 
A. Background 

Vertically oscillated granular layers provide an impor- 
tant testbed for granular research. When shaken verti- 
cally, fiat layers of grains exhibit convection [l|, cluster- 
ing [ah shocks [H, steady-state flow fields far from the 
platejH, floating particle clusters @, and standing- wave 
pattern formation Q. 

In this paper, we model granular media by numerically 
solving a set of time-dependent continuum equations for 
the rapid flow of dissipative particles in three-dimensions. 
We use these simulations to investigate the relationship 
between shocks and standing-wave patterns in vertically 
oscillated granular layers. 



B. Granular hydrodynamics 

A successful theory of granular hydrodynamics would 
allow scientists and engineers to apply the powerful meth- 
ods of fluid dynamics to granular flow. Despite ex- 
perimental 0-Q and computational [1, [Iol - [l3j evidence 
demonstrating the potential utility of hydrodynamics 
models for grains, a general set of hydrodynamic gov- 
erning equations is not yet recognized for granular media 

Several proposed rapid granular flow models use equa- 
tions of motion for continuum fields: number density n, 
velocity u, and granular temperature T (§T is the aver- 
age kinetic energy due to random particle motion) fl8l — 
|20| . In one approach, particle interactions are modeled 
with binary, inelastic hard-sphere collision operators in 
kinetic theory to derive continuum equations to Euler 
[2lj . Navier-Stokes H3, and Burnett 23} order. We use 
three-dimensional (3D) simulations of continuum equa- 
tions to Navier-Stokes order to investigate shocks and 
standing-wave patterns in shaken granular layers. 



C. Standing wave patterns in oscillated granular 
layers 

A granular layer of depth H atop a plate that is os- 
cillated sinusoidally in the direction of gravity with fre- 
quency / and amplitude A will leave the plate at some 
time during the oscillation cycle if the maximum accel- 
eration of the plate a max = A (2-rrf) 2 is greater than the 
acceleration of gravity g. In other words, the layer leaves 
the plate if the nondimensional accelerational amplitude 
T = a m ax/ g exceeds unity. 

After leaving the plate, the layer dilates above the 
plate, and then compresses when it collides with the 
plate later in the cycle. When the dimensionless ac- 
celerational amplitude T is larger than a critical value 
Tc, standing- wave patterns spontaneously form in the 
layer. These patterns are subharmonic with respect to 
the plate, repeating every 2// [6]. Depending on the 
nondimensional accelerational amplitude T and dimen- 
sionless frequency /* = fy/H/g, various subharmonic 
standing waves have been found, including stripe, square, 
and hexagonal standing- wave patterns [6j . 



D. Shocks in oscillated granular layers 

If the Mach number Ma (the ratio of the local mean 
fluid speed to the local speed of sound) is greater than 
unity at the point where a fluid encounters an obsta- 
cle, a compression wave front is formed near the object 
and steepens to form a shock. A distinguishing feature of 
granular materials is that granular flows reach supersonic 
speeds under common laboratory conditions [3 H, fl8j - 
Therefore, although shocks are formed in ordinary gases 
only under extreme conditions, shock formation is com- 
monplace in granular media. 

Experiments Q and simulations [ll|, HiHIsI demon- 
strate that shock waves form in shaken granular layers as 
the layer contacts the plate. Previous investigations have 
generally either focused on shock propagation [H, HU 
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or pattern formation [6j, |l2|, ll3| or considered them as 
separate phenomena coexisting in shaken layers [TTj . In 
this paper, we use continuum simulations to investigate 
the interactions between shocks and the standing-wave 
patterns formed in this system. 



E. Model system 

We simulate a layer of grains on an impenetrable plate 
which oscillates sinusoidally in the direction of gravity. 
The layer depth at rest is approximately H = 5.4er, where 
the grains are modeled as identical, frictionless spheres 
with diameter a, mass M, and coefficient of restitution 
e = 0.7. In this paper, we study patterns and shocks 
as a function of nondimensional frequency /*, while the 
dimensionless accelerational amplitude T = 2.20 is held 
constant. Previous simulations [H, [l]| have shown that 
this value exceeds the critical accelerational amplitude 
Tc indicating that we should expect to see standing- wave 
patterns for a variety of frequencies. 

and molecular dynamics 
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Previous experiments 
(MD) simulations [28[ have shown that friction between 
grains plays a role in these patterns. Experimentally, 
adding graphite to reduce friction decreased Tc and pre- 
vented the formation of stable square or hexagonal pat- 
terns found for certain ranges of /* and T in experi- 
ments without graphite [27j . Similarly, MD simulations 
with friction between particles have quantitatively repro- 
duced stripe, square, and hexa gon al subharmonic stand- 
ing waves seen experimentally [29(, but MD simulations 
without friction yield only stable stripe patterns and dis- 
play a lower Tc [Hj]. In this study, we investigate stripe 
patterns in continuum simulations of frictionless parti- 
cles. To investigate other patterns such as squares or 
hexagons, simulations would have to include friction be- 
tween particles. 

Experiments [3(| and simulations [lU, [l3|, [3l[ indicate 
that fluctuations due to individual grain movement play 
a larger role in granular media than do thermal fluc- 
tuations in ordinary fluids. Fluctuating hydrodynam- 
ics (FHD) theory models these fluctuations by adding 
noise terms to the Navier-Stokes equations |32h34| . In 
previous simulations with FHD noise terms, the critical 
value of pattern onset Tc was consistent with molecu- 
lar dynamics simulations, while continuum simulations 
without these fluctuations exhibited pattern onset at Tc 
approximately 10% lower than that found in MD simula- 
tions ^2l.[l3j. Above Tc, however, simulations both with 
and without these fluctuations exhibited standing-wave 
patterns with wavelengths consistent with a dispersion 
relation [35j | found experimentally for a range of shaking 
frequencies. In this paper, we do not include FHD noise 
terms; we investigate patterns for accelerational ampli- 
tude greater than Tc determined from simulations with 
and without FHD. 

We use continuum simulations to investigate the dy- 
namics of this system including pattern formation and 



shock propagation. Section |TT] describes the methods we 
use to simulate and analyze oscillated layers. Sec. IIIII 
examines the dynamics of shocks and standing- wave pat- 
terns formed in this system at a fixed dimensionless os- 
cillation frequency /* = 0.25, and Sec. IIVI examines how 
these patterns and shocks change when this frequency is 
varied. We present our conclusions in Sec. [V] 



II. METHODS 

A. Continuum equations 

We use a continuum simulation previously used to 
model shock waves in a granular shaker 26]. Our sim- 
ulation numerically integrates continuum equations of 
Navier-Stokes order proposed by Jenkins and Richman 
22] for a dense gas composed of frictionless (smooth), 
inelastic hard spheres. 

This model yields hydrodynamic equations for number 
density n (or equivalently, volume fraction v — ^nc 3 ), 
velocity u, and granular temperature T: 



du 

~dt 



u • Vu I = V • P — ngz, 



(1) 



(2) 



|n^+u-Vr) =-V-q + P:E- 7 , (3) 

where the components of the symmetrized velocity gra- 
dient tensor E are given by: = 5 (djUi + diUj) . The 
components of the stress tensor P are given by the con- 
stitutive relation 



P — 



-p+ (A - -fi)E kk 



Sij+2fiEij, (4) 



and the heat flux is calculated from Fourier's law: 

q = -kVT. (5) 

To calculate the pressure, we use the equation of state 
and radial distribution function at contact proposed by 
Goldshtein et al. [3(| to include both dense gas and in- 
elastic effects: 



p = nT [1 + 2(1 + e)G{v)\ , 

where 

G(v) = vg (v), 
and the radial distribution function at contact, go, is 



9o(v) 
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(6) 
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where u max — 0.65 is the 3D random close-packed volume 
fraction. 

These equations differ from those for a compressible, 
dense gas of elastic particles by the energy loss term 



7 = il(l-e 2 ) 



nT 3 / 2 



G{y), 



(9) 



which arises from the inelasticity of collisions between 
particles. The bulk viscosity is given by 
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the shear viscosity by 
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and the thermal conductivity by 
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Other hydrodynamic models include modihcations 
such as more accurate expressions for kinetic coeffi- 
cients which incorporate high density corrections (37l . l38j 
and equations of state which allow for volume fractions 
gre ater than the random-close-packed limit used here 
[391 ] . However, the equations shown here have previously 
been used to separately examine shocks [26| and patterns 
[l2l ll3T | in shaken layers and have demonstrated quantita- 
tive agreement with experiments and molecular dynam- 
ics simulations. Therefore, we use these equations in this 
study to investigate the relationship between shocks and 
pattern formation; in principle, the equations could be 
modified to implement other forms of the constitutive 
relations. 



B. Simulation method 

We integrate these hydrodynamic equations to find 
number density, momentum, and granular temperature, 
using a second order finite difference scheme on a uniform 
grid in 3D with first order adaptive time stepping [26| . 
In these simulations, the granular fluid is contained be- 
tween two impenetrable horizontal plates at the top and 
bottom of the container, where the lower plate oscillates 
sinusoidally between height z — and z = 2A and the 
ceiling is located at a height L z above the lower plate. Pe- 
riodic boundary conditions are used in the horizontal di- 
rections x and y to eliminate sidewall effects. Simulations 
were conducted in a box of size L x = 168cr, L y = 10a, 
and L z = 160o\ This orientation causes stripes to form 
parallel to the j/-axis. The numerical methods, boundary 
conditions at the top and bottom plate, and grid spa cing 
are the same as used in previous studies of shocks [26| 
and patterns [IH, [l3| . 



III. DYNAMICS OF SHOCKS AND STANDING 
WAVES 

A. Standing wave patterns 

Experimental investigations of shaken granular layers 
have shown that above a critical acceleration of the plate 
r<7, standing wave patterns form spontaneously. These 
patterns oscillate subharmonically, repeating every 2/f, 
so that the location of a peak of the pattern becomes a 
valley after one cycle of the plate, and vice versa @. 

Continuum simulations produce standing wave pat- 
terns for T = 2.2 and /* = 0.25 (Fig. [T). Beginning 
with a flat layer above the plate with small amplitude 
random fluctuations, the simulation ran for 250 cycles of 
the plate until the layer reached a periodic state. Snap- 
shots from various times during the next two cycles of the 
plate are shown in Fig. [T] Alternating peaks and valleys 
form a stripe pattern which oscillates at f/2 with respect 
to the plate oscillation; a location in the cell which rep- 
resents a peak during one cycle will become a valley the 
next cycle, and then return to a peak on the following 
cycle. We examine the wave patterns at various times in 
the cycle ft for two cycles of the plate. 



1. ft — 0: The layer collides with the plate. 

At ft — the container is at its minimum height. The 
bulk of the layer is dilute and four clear peaks and four 
valleys are evident in Fig. [T] (note that one of the peaks 
wraps around the edges of the cell since the boundary 
conditions are periodic). The material in each peak is 
falling towards the plate, although the bottom of the 
layer has begun to contact the plate and is beginning 
to form a pile on the plate. The thin layer colliding with 
the plate is characterized by a higher volume fraction v 
than the material in the peaks. 



2. ft = 0.25: The layer piles on the plate. 

At ft = 0.25, the plate has begun move upward as the 
layer falls to meet it. At this time, most of the layer is 
piling on the plate, forming a high volume fraction region 
near the plate. Remnants of the peaks and valleys from 
the previous cycle are still visible, although the layer is 
flattening on the plate, with the height difference between 
peaks and valleys much smaller than at ft = (Fig. [T]). 



3. ft — 0.50: The layer begins to leave the plate. 

As the container reaches its maximum downward ac- 
celeration at the top of its oscillation, the layer begins 
to leave the plate. Although the bottom of the layer 
is still in contact with the plate at this time, the layer 
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is expanding, as evidenced by the fact that it is visibly 
more dilute in this picture than in the previous snapshot 
(Fig. [I]). As the layer leaves the plate, new peaks and 
valleys develop, with peaks forming in the valleys from 
the previous cycle, and vice versa. 
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FIG. 1: A side view of a layer of grains, showing volume 
fraction v at a slice y — 5er parallel to the x—z plane at various 
times ft during two cycles of the plate. Empty space (u — 
0) is white, random close-packed volume fraction (c ma x = 
0.65) is black, and the color increases from white to black 
through various shades of gray as volume fraction increases. 
The plate is represented as a thick, horizontal black line. The 
cell extends to a height of 160a above the lower plate, but 
the figures show only z < 50<r, since the density is quite low 
above this height. 



4- ft = 0.75: The layer is off the plate. 

At ft = 0.75, the plate is moving downward leaving 
a gap between the layer and the plate. The mass of the 
layer has almost entirely left the plate, and the layer has 
expanded. The new peaks and valleys have become quite 
distinct, with a large portion of the layer in the peaks 
and very little material in the valleys connecting them 
(Fig. □]). 



5- 1.0 < ft< 2.0: The cycle repeats. 

At ft — 1.0, the plate has undergone one full oscillation 
and has returned to its lowest point in the cycle. The 
features of this point in the cycle are similar to those 
at ft = 0, except the peaks and valleys have reversed 
location. Aside from the peaks and valleys reversing, the 
next cycle exhibits the same features as the previous cycle 
at all times during the cycle. By ft — 1.75 the peaks and 
valleys are back to their original location, demonstrating 
the subharmonic nature of these patterns. 



B. Shocks 

Previous experiments [|[ and simulations [IJ HHIll 
indicate that shocks are formed in vertically oscillated 
granular layers during each collision of the layer with 
the plate. A prerequisite for shock formation is that the 
local Mach number of the flow be greater than unity with 
respect to the object causing the disturbance. Therefore, 
we calculate the speed of sound using a relation derived 
from the equation of state Eq. © [40| : 



1 + t;X + -rtfl - 

6 cm ov 



ITx 



(13) 



where \ — 1 + 2(1 + e)G{v). The Mach number Ma is 
calculated as the local flow speed divided by the local 
speed of sound Ma — |u| /c. 

From Fig. [TJ we see that at ft = 1.00, a portion of 
the layer has begun to compress on the plate, while some 
of the material is still falling towards the plate. This is 
therefore a time at which we would expect to see shock 
formation. To examine the shock in detail, we plot one- 
dimensional profiles of volume fraction v and Mach num- 
ber Ma (here v is scaled by a factor of 10 to fit on the 
same scale as Ma) as functions of height z from the plate 
in Fig. H 
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FIG. 2: Mach number Ma (dashed line) and rescaled volume fraction 10^ (solid line) as functions of height z (ordinate) at two 
times ft in the oscillation cycle. In each case, the plate is shown as a horizontal solid black line. The left column shows Ma 
and lQv at the horizontal location y — 5a, x — 16a at times (a) ft = 1.0 and (c) ft — 1.1 in the plate's oscillation, while the 
right column shows Ma and lOv at the horizontal location y = 5a, x — 32a at times (b) ft — 1.0 and (d) ft — 1.1. A shock is 
present in each case, as indicated by a discontinuity in the volume fraction derivative and a sharp increase in Ma moving away 
from the plate. 



At ft = 1.00, the first peak of the pattern occurs at 
the horizontal location 16a < x < 20a (c/Fig. [1]). We 
choose two horizontal locations to examine; the location 
y = 5a, and x — 16a is shown in the left column of 
Fig. [51 corresponding to the first peak of the pattern; 
the location y = 5a, and x = 32a is shown in the right 
column, corresponding to a location away from the peak. 

At the horizontal location y — 5a, x = 16a correspond- 
ing to a peak of the wave pattern, the volume fraction 
approaches the close-packed value v ss 0.6 near the plate 
at ft = 1.00 (Fig. [5^). As height increases, v smoothly 
decreases for < z < la. The Mach number with respect 
to the plate Ma is low throughout this region, as this part 
of the layer is compressed on the plate and thus matches 
the velocity of the plate at height z = 0. However, at 
z ~ la, there is a discontinuity in the derivative of v, 
and a sharp increase in Mach number such that Ma « 3 
for z > 8a. This sharp increase in Mach number from 
subsonic to supersonic, corresponding to a discontinuity 
in the derivative of volume fraction at the same height, 
indicates that this is the location of a strong shock front. 

While the above description applies to x — 16a (the 
horizontal location of the first peak), a shock also forms 
in the valleys. This can be seen in Fig. ^jp, which shows 
lOi' and Ma at the same time, but at horizontal location 



x = 32a. Here we again see a high density, subsonic 
region near the plate, with a shock separating this region 
from a lower density, supersonic undisturbed region that 
is still falling towards the plate. However, the volume 
fraction near the plate is smaller than it was in Fig. [5^, 
and the shock profile is quite different. Also, the shock 
front at x = 32a is closer to the plate than it is at x = 
16a. 



By ft = 1.10, the shock has moved away from the plate 
and it continues to develop as it moves through different 
parts of the layer (Fig. The shock front in the valleys 
is also developing and moving away from the plate (see 
Fig. [2ji) , although the location and profile of the shock 
is different at these different horizontal locations. 



Thus, shocks form with each collision of the layer with 
the plate, and the shock profile changes and develops as 
the shock moves through the layer. Shocks form at each 
horizontal location in the cell; however the shock front is 
non-uniform horizontally due to the differences in layer 
depth between the peaks and the valleys of the pattern. 
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C. Interaction between shocks and patterns 

There is therefore a relationship between the shocks 
and the standing-wave patterns formed in the layer. If 
the layer were perfectly flat, the shock front would be uni- 
form horizontally. However, as there is horizontal varia- 
tion in layer depth, a non-uniform shock front forms. 

While the system is driven by vertical motion of the 
plate, the patterns are characterized by horizontal varia- 
tion between peaks and valleys. Periodic horizontal mo- 
tion is required to produce the subharmonic oscillation. 
This horizontal flow can be thought of as a "sloshing" 
motion, in which material from the layer flows out of 
some regions and into others, creating peaks and valleys. 
With the next collision of the layer with the plate, the 
flow direction reverses, creating peaks where there were 
valleys, and vice versa. In this section, we show that 
strong pressure gradients created by the shock drive the 
sloshing motion of the layer. 

To study this interaction, we investigate the layer prop- 
erties near the time of shock formation. According to 
Fig.[T] most of the layer is off the plate at ft = 0.75, while 
the layer is compressed on the plate and the peaks and 
valleys have flattened significantly by ft = 1.25. In Fig. [3J 
we show snapshots of the the volume fraction v (left col- 
umn) and the dimensionless pressure p (right col- 
umn) at five different times 0.8 < ft < 1.2. We display 
side views of the layer at y — 5a in (Fig [3]), showing only 
the horizontal range < x < 84a and the vertical range 
< z < 25cr to better view the interaction between the 
layer and the plate. This range shows two wavelengths of 
the standing wave pattern and corresponds to the bottom 
left quadrant of the range plotted in Fig. [TJ 

In order to examine the sloshing motion of the layer, we 
also show the dimensionless x-component of the pressure 

gradient || {^fj^j (left column) and the dimensionless 

horizontal velocity u x ^fga (right column) at the these 
same times and locations in Fig. |H We now proceed 
to examine the horizontal variation in flow at these five 
times as the layer collides with the plate. 

1. ft~ 0.8 : The layer is off the plate. 

At ft = 0.8, the layer is mostly off the plate, and peaks 
and valleys are clearly visible in the left column of Fig. |3j 
In the right column of Fig. |3j we show the dimensionless 
pressure p (jfjj) ■ At this time, the layer is more dilute 
than it will be at later times in the cycle, and the pres- 
sure is uniformly low throughout the layer. There are no 
strong pressure gradients visible, as can be seen in the 
left column of Fig. 21 where we show the dimensionless 

x-component of the pressure gradient |^ (ivfg) • 

In the right column of Fig. SI we can see that u x > 
(the flow is moving towards the right) for < x % 16a\ 
there is a small range of nearly zero horizontal flow rate 



near 16a < x < 20cx; then u x > (flow is moving towards 
the left) for 16a < x < 34c Since the peak of the pattern 
is located at 16a < x < 20er at this time (c/Fig.[3j, grains 
are flowing towards the peak from both the left and the 
right. Therefore, the flow is still removing material from 
the valleys and moving it towards the peak; the difference 
between peaks and valleys is growing at this time. 

2. ft ~ 0.9 : The layer begins to contact the plate. 

At ft = 0.9, the bottom of the layer has just begun 
to contact the plate at the horizontal location of the two 
peaks, but the material between the peaks has not yet 
contacted the plate (Fig. [3J . The layer looks similar to 
its appearance at ft = 0.8, except that in a very small 
region near the plate, the pressure has begun to slightly 
increase where the layer is starting to contact the plate. 

This pressure increase at the two points of contact 
yields a horizontal pressure gradient with maximum pres- 
sure near the peak location. For instance, the reddish 
color found near the plate for 12ct < x < 18cr in the 
left column of Fig. [4] at this time indicates that the pres- 
sure is increasing from left to right; the bluish color for 
18a < x < 26a in the same picture indicates that the 
pressure decreases from left to right. 

The picture in the right column for this time shows 
that the flow is still moving towards the peak through 
the bulk of the layer at this time. However, the pressure 
gradient near the plate tends to drive flow in the oppo- 
site direction: from high pressure (near the peak) into 
low pressure (near the valley). Therefore, the horizontal 
flow velocity near the plate has begun to slow (leaving a 
white space near the plate where the pressure gradient 
is strongest) and is even beginning to reverse direction 
(note the light blue near the plate at the location of the 
second peak). 

3. ft » 1.0 : A shock forms. 

By ft — 1.0, more of the material is colliding with the 
plate, and the pressure has noticeably increased at the 
two main contact points between the layer and the plate 
(Fig. [3J . This increase in layer density and pressure near 
the plate marks the formation of a shock (cf Fig. 

As discussed in Sec. lIIIBl these shocks are not uniform 
horizontally; larger layer depth near the peaks leads to 
larger pressure at these locations. Thus, shock formation 
leads to pressure maxima near the plate at the horizontal 
location of the pattern peaks with strong pressure gradi- 
ents as seen in the left column of Fig. @] 

These pressure gradients drive the flow in the direc- 
tion opposite the pressure gradient; thus the direction of 
flow near the plate at this time is clearly reversing (see 
the right column of Fig. 2]). For much of the layer, the 
flow is still away from the valleys and towards the peaks. 
Near the plate in the region of strong pressure gradients, 
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FIG. 3: A side view of a layer of grains, showing volume fraction v (left column) and dimensionless pressure p (right 

column), at a slice y — 5a parallel to the x — z plane at various times ft. Although the full cell has a length 168a in the 
x-direction and a height 160a in the ^-direction, this figure only displays the horizontal range < x < 84a and the vertical range 
< z < 25a to show a closer view of the collision of the shock with the plate. In the left column, empty space (y = 0) is white, 
random close-packed volume fraction [y = 0.65) is black, and the color increases from white to black through shades of gray 
as volume fraction increases. The plate is represented as a thick, horizontal black line. In the right column, low dimensionless 
pressure is white, high dimensionless pressure is black, and pressure increases from white to black through shades of gray. 



however, the flow has reversed itself and is now flowing 
away from the peaks and towards the valleys. This sep- 
aration of the flow into two distinct domains represents 
the formation of a shock near the plate. 

4- /t ~ 1.1 : The shock propagates. 

This shock then travels upward through the layer, away 
from the plate. By ft = 1.1, this shock is well-developed 
and has propagated through most of the layer, separating 
a high- density and pressure region near the plate from 
a low- density and pressure region above the shock front 
(Fig.©. 

The pressure gradient is quite strong near the plate 



(Fig. 2]) and most of the layer is behind the shock. The 
flow behind the shock has reversed direction, while the 
material ahead of the shock has not yet changed direc- 
tion. There is still some material in the peaks that is still 
falling towards the plate ahead of the shock (c/Fig. [3]), 
but its horizontal velocity is small, and the pressure is 
relatively low in this region (see Fig. @|. 

5. ft fa 1.2 : The flow has reversed. 

By ft — 1.2, nearly the entire layer is behind the shock 
and the layer has noticeably flattened (Fig. [3]), although 
peaks and valleys are still visible. As the shock prop- 
agates into the low-density region above the layer, the 
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FIG. 4: A side view of a layer of grains, showing dimensionless horizontal pressure gradient ("Pg) (^ e ^ column) and 
dimensionless horizontal velocity u xy fga (right column) at a slice y = 5a parallel to the x-z plane at various times ft. 
Although the full cell has a length 168<r in the s-direction and a height 160cr in the z-direction, this figure only displays the 
horizontal range < x < 84a and the vertical range < z < 25a to show a closer view of the collision of the shock with the 
plate. In the left column, pressure decreasing left to right (J^ < 0) is shown in shades of blue, ^ = is white, and pressure 
increasing left to right (|£ > 0) is shown in shades of red. In the right column, flow towards the left (u x < 0) is shown in 



shades of blue, zero horizontal velocity (y = 0) is white, and flow towards the right (u x > 0) is shown in shades of red. In both 
columns, the plate is shown as a thick, horizontal black line. 



layer cools behind the shock due to inelastic collisions. 
This causes the pressure to decrease slightly at ft = 1.2 
as compared to ft = 1.1. 

The pressure gradient is still driving the flow away from 
the peaks, although this gradient is not as strong as in 
at ft = 1.1 (Fig. At this time, nearly the entire layer 
has reversed direction and is now flowing away from the 
peaks and towards the valleys. 

As this flow continues, the layer becomes flat, and then 
develops peaks where there were previously valleys, and 
vice versa. Thus, the horizontal pressure gradients cre- 
ated by the nonuniform shock front drive the flow to re- 
verse itself. It is this reversal that creates the sloshing 
motion visible in the subharmonic oscillation. 



IV. FREQUENCY DEPENDENCE OF SHOCKS 
AND STANDING WAVES 



In Sec. IIII1 we studied the interaction between shocks 
and patterns in layers oscillated with a particular non- 
dimensional frequency /* = 0.25. In this section, we 
investigate frequency dependence by holding dimension- 
less accelerational amplitude V = 2.2 constant, while 
varying dimensionless frequency /*. Experiments have 
shown that wavelength A depends on the frequency of 
oscillation [H, |4l|. For a range of layer depths and 
oscillation frequencies, experimental data for frictional 
particles near pattern onset were fit by the function 
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A* = 1.0 + l.l/*-l-32±0.03 5 where y = X j H ^ We 

examine the correlation between changes in shock proper- 
ties with changes in pattern wavelength throughout this 
frequency range. For each frequency, we start from an 
artificial flat layer, then simulate 250 cycles of the plate 
to allow the layer to reach an oscillatory state. Data is 
taken from the next six cycles of the plate. 



well as the time during the cycle at which the collision 
takes place, higher frequencies correspond to lower Ma 
at collision throughout this range. 



B. Higher Mach number at collision produces 
stronger shocks 



A. Maximum Mach number varies inversely with 
driving frequency 

Shock properties depend on the Mach number of the 
layer with respect to the plate during collision. In Fig. [5j 
we plot the maximum Mach number Ma mas of the layer 
with respect to the plate as a function of dimensionless 
driving frequency /*. To calculate Ma max , we find the 
Mach number Ma (x, t) at each location x in the cell 
during an oscillation cycle. To ensure that we are looking 
at the Mach number of material in the layer, rather than 
in the low-density region above the layer, at each time 
t, we find the highest Mach number corresponding to at 
least 1% of the total mass of the layer. We then define 
Afa malt as the maximum Mach number found at any time 
during a cycle. 



30 



Ma 
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Pressure decreases rapidly when moving from the re- 
gion behind a shock, across the shock front, to the undis- 
turbed region ahead of the shock. We examine this 
pressure change by calculating the magnitude of the z- 



component of the pressure gradient 



dp 



at each 

(*) 



time t and location in the cell x. We then find 

as the highest value corresponding to at least 1% of the 
total mass of the layer. Finally, we calculate the time 



average of these layers 



Oz 



, nondimensionalize by 



a factor of jj^, and plot this as a function of Ma max 
in Fig. [5a] As the maximum Mach number Ma max in- 
creases, the time-averaged maximum vertical pressure 
gradient increases monotonically. In other words, the 
higher the relative Mach number of the layer with re- 
spect to the plate, the sharper the pressure drop across 
the shock. Note that hydrostatic pressure would cause 
pressure to vary with depth even in static layers, consis- 
tent with the non-zero intercept of Fig. I6al 

As discussed in Sec. Mil although the shock is produced 
by collision with a vertically oscillating plate, horizontal 
pressure variation develops as a result of the peaks and 
valleys in the layer. Shocks with stronger vertical pres- 
sure gradients correspond to stronger horizontal pressure 
gradients as well (Fig. I6bl) . A flat layer with no patterns 
would be expected to have pressure variation with depth 
but not with horizontal location. This is consistent with 
the apparent non-zero intercept in Fig. I6bl 
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FIG. 5: The maximum Mach number Ma mB of the layer with 
respect to the plate found at any time during an oscillation 
cycle as a function of the dimensionless driving frequency /*. 
We calculate Ma max for each of six oscillation cycles; points 
are the average of the six cycles and error bars represent the 
standard deviation. 

Note that Ma max monotonically decreases as /* in- 
creases (Fig. [5J. For fixed T = 2.20, lower oscillation 
frequency corresponds to a higher maximum plate veloc- 
ity and higher oscillation amplitude. Layers with lower 
/* will generally reach a higher maximum height when 
they leave the plate, and will have a larger speed relative 
to the plate when they collide later in the cycle. Al- 
though the relative Mach number between the layer and 
the plate depends on the speed of sound in the layer as 



C. Horizontal pressure gradients drive horizontal 
velocity 

As strong pressure gradients develop in the layer, there 
will be a tendency for the material to flow from high 
pressure to low pressure; i.e. the direction of flow ve- 
locity will be opposite the pressure gradient (c/Fig. 0}. 
We calculate the average flow speed in the cc-direction at 
all times and all locations in the cycle Figure [7] 

shows that the nondimensionalized average flow speed 
in the x— direction (|u x |) / ' \fg& increases monotonically 
as the maximum dimensionless x-component of the pres- 
sure gradient found anywhere in the cell averaged over 



all times in the cycle 



Mg 



increases. Thus, the 



stronger the maximum horizontal pressure gradients pro- 
duced across the shock are, the faster the average hori- 
zontal flow speed becomes. 
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max 





FIG. 7: The nondimensional average flow speed in the x- 
direction (ImxI) / as a function of the dimensionless max- 
imum ^-component of the pressure gradient found anywhere 
in the cell averaged over all times in the cycle (| ff | max ) jj^- 
Points are calculated as the average of six oscillation cycles of 
the plate, and error bars represent the standard deviation of 
these six cycles. 



wavelength yields inherent uncertainty in the wavelength 
that would be selected in an infinite box. 

Faster flow in the x-direction corresponds to more hor- 
izontal motion of the particles throughout an oscillation 
cycle. Therefore, the wavelength (and therefore the dis- 
tance between peaks) increases monotonically as the av- 
erage horizontal flow speed increases, (see Fig. [3]). 



FIG. 6: The dimensionless maximum ^-component (a) of the 
pressure gradient found anywhere in the cell averaged over all 
times in the cycle ( I i 2 1 ) -S- as a function of the maxi- 

J \ I oz I max / Mg 

mum Mach number Ma max of the layer with respect to the 
plate found at any time during an oscillation cycle, and the 
maximum ^-component (b) of the pressure gradient found 
anywhere in the cell averaged over all times in the cycle 

(Iff I max) Tig as a f unc tion of the maximum ^-component 
of the pressure gradient found anywhere in the cell averaged 
over all times in the cycle ) 4i- ■ We calculate values 

J \ I az I max/ Mg 

for each of six oscillation cycles; points are the average of the 
six cycles and error bars represent the standard deviation. 
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D. Horizontal velocity produces standing waves 

For the dimensionless frequency /* = 0.25 examined in 
Sec. IIII1 four wavelengths fit in a box of size 164cr in the 
x-direction, yielding a wavelength of 41a (Fig. [I}. We 
calculate the dimensionless wavelength A* = X/H, and 
plot it as a function of nondimensional average flow speed 
in the ^-direction |) / ^fgd for each frequency in our 
range (Fig. [8j). Due to the periodic boundary conditions 
and finite size of the box, an integer number of waves 
must fit in the box. This finite size effect of quantized 



1 (\u x \)/y/ga 2 

FIG. 8: Nondimensional wavelength A* = X/H as a func- 
tion of nondimensional average flow speed in the x-direction 
(\ux\) /yfgv- Points are calculated as the average of six os- 
cillation cycles of the plate. The error bars for A* = X/H 
are calculated exclusively from discretization due to periodic 
boundary conditions in a finite-size box. The error bars for 
(\u x \) / \fgo represent the standard deviation of these six cy- 
cles. 
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E. Dispersion relation 

Summarizing the results from this section thus far, we 
find that higher frequencies of oscillation yield smaller 
Mach numbers between the layer and the plate during 
collision (Fig. [5]). Larger values of the Mach number 
produce stronger shocks, as reflected in larger vertical 
and horizontal changes in pressure across the shock front 
(Fig. Larger pressure gradients produce faster hori- 
zontal flow (Fig. [7]) which in turn corresponds to larger 
horizontal distances between peaks of the pattern, and 
thus a larger wavelength (Fig. HJ). 

Given this chain of reasoning, it follows that as the 
Mach number of the layer with respect to the plate in- 
creases, wavelength of the pattern increases monotoni- 
cally throughout the range. As Mach number increases 
with decreasing shaking frequency, wavelength should 
decrease monotonically as the oscillation frequency in- 
creases. Figure M shows that this is indeed the case 
throughout the range of frequencies we have simulated. 
Additionally, the wavelengths found in our simulation 
throughout this range are consistent with the dispersion 
relation previously found by fit to experimental data [35j j . 




.1 .2 .3 .4 f * 

FIG. 9: Dimensionless wavelength A* as a function of dimen- 
sionless frequency /* . Data from our simulations are shown as 
points, with error bars calculated exclusively from discretiza- 
tion due to periodic boundary conditions in a finite-size box. 
The dominant wavelength found in our simulations fits quite 
well the dispersion relation A* = 1.0 + ij/*- 1 - 32 ± - 03 (solid 
line) found as a fit to previous experimental data [35| ■ 
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V. CONCLUSIONS 

We have simulated vertically oscillated layers of granu- 
lar media by numerically solving a proposed set of granu- 
lar hydrodynamic equations to Navier-Stokes order. We 
have shown that these continuum simulations can de- 
scribe important aspects of shocks and patterns in gran- 
ular materials. In our simulations, layers spontaneously 
form subharmonic standing-wave patterns in which the 
wavelength A depends on the frequency / of the oscil- 
lation of the plate. The wavelengths of the patterns 
formed in our simulation agree with the dispersion re- 
lations found in previous experiments. 

As these standing waves oscillate subharmonically, 
shocks are formed with each collision of the layer with 
the plate. We have analyzed these shocks and found that 
variation in the layer depth produces a shock front that 
is not uniform horizontally. This horizontal variation in 
the shock front produces horizontal components of flow 
velocity, causing grains to move from high pressure to low 
pressure. Therefore, despite the fact that the plate oscil- 
lates vertically, the shocks produced during collision be- 
tween the layer and the plate drive the horizontal sloshing 
motion that characterizes the standing-wave patterns. 

We have simulated shocks in layers oscillated at vari- 
ous frequencies, and have shown that the strength of the 
shock varies depending on the oscillation frequency. We 
have also established that there is a relationship between 
pattern wavelength and the Mach number of the layer 
with respect to the plate at the time of collision. Specifi- 
cally, the horizontal and vertical components of the pres- 
sure gradient both increase monotonically as the Mach 
number increases. The horizontal flow speed increases 
monotonically with the strength of these pressure gradi- 
ents, leading to longer wavelengths. 

Therefore, throughout the range we studied, the wave- 
length of the pattern increases monotonically with the 
maximum Ma of the layer with respect to the plate. This 
analysis indicates that shocks play a significant role in the 
dynamics of the standing-wave patterns formed in these 
oscillating layers. 
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